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For a Bose-Einstein condensate placed in a rotating trap and confined in the z axis, we set a 
framework of study for tlie Gross-Pitaevskii energy in tfie Tliomas Fermi regime. We investigate an 
asymptotic development of the energy, the critical velocities of nucleation of vortices with respect 
to a small parameter e and the location of vortices. The limit e going to zero corresponds to the 
Thomas Fermi regime. The non-dimensionalized energy is similar to the Ginzburg-Landau energy for 
superconductors in the high-kappa high-field limit and our estimates rely on techniques developed for 
this latter problem. We also take the advantage of this similarity to develop a numerical algorithm 
for computing the Bose-Einstein vortices. Numerical results and energy diagrams are presented. 
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I. INTRODUCTION 

Since the first experimental achievement of Bose- 
Einstein condensates in atomic gases in 1995, many prop- 
erties of these systems have been studied experimentally 
and theoretically, and particularly the existence of vor- 
tices H, H |, |, |, |, 0, |]. A way to create vortices 
consists in rotating the trap confining the atoms: for suf- 
ficiently large velocities, it becomes energetically favor- 
able to have vortices in the system. Theoretical stud- 
ies of this type of experiments have often been made 
in the framework of the nonlinear Schrodinger equation 
or Gross-Pitaevskii equation, well known for superflu- 
ids, but which provides a very good description of Bose- 
Einstein condensates: it is assumed that the N particles 
of the gas are condensed in the same state for which the 
wave function (j) minimizes the Gross-Pitaevskii energy. 
By introducing a rotating frame at the angular velocity 
Q, = Sle^, the trapping potential becomes time indepen- 
dent, and the wave function minimizes the energy 

^2 



N 



g3D|0|"-?ij^-(i0,V0xx), (1) 



under the constraint / If/ip = 1. Here, for any complex 
quantities u and v and their complex conjugates u and v, 
(u, v) = {uv + uv)/2. The terms in the energy correspond 
to the kinetic energy, the trapping potential energy, the 
interaction energy and the inertial due to the change of 
frame. 
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When the atoms are strongly confined in the z— axis, 
the situation can be simplified to a two dimensional prob- 
lem where the wave function -i/; depends only on x = (x, y) 
and it minimizes 

a—x,y 

N 

+-g\^Pf-hn-it^,V^Px^), (2) 

where g — g^Dimuj z / 27:1%) . The constraint J 1-01^ = 1 is 
also imposed. Our study was at first motivated by the 
work of Castin and Dum j2| who have studied the equi- 
librium configurations by looking for the minimizers in a 
reduced class of functions for the 2D case and have done 
numerical computations in 2 and 3D. Their analysis is in 
the Thomas-Fermi regime, where the mean interaction 
energy per particle is larger than fiujx.y- 

Our aim is to provide a mathematical framework for 
a rigorous study of the energy 820 and its minimizers in 
the Thomas Fermi limit. We first observe that this en- 
ergy has a striking similarity with the high-kappa, high- 
field limit of the Ginzburg-Landau free energy used in 
the modeling of superconductors. Thus, we expect the 
energy £20 will develop similar behavior as those for the 
Ginzburg-Landau energy studied in [|ll|, In par- 

ticular, the results obtained in the context of Ginzburg- 
Landau energy may be applied to £20 , in the Thomas 
Fermi regime, to yield an asymptotic development of the 
energy as well as the critical velocities for the nucleation 
of vortices and the location of these vortices. Due to the 
close resemblance, we will not be concerned with the de- 
tailed derivations in this paper, but rather focus on the 
conclusions one can draw from the asymptotic develop- 
ments. To our knowledge, some of our estimates to be 
presented later have not been given in the literature pre- 
viously. Let us point out that our method to compute 
the energy is also new in this context and very different 
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from the ones in, for instance, |g, ^. 

We define the characteristic length d = {fi/mujxY/'^ 
and set 



Note that critical point u of is a solution of 
- Au - 2i{n X r).Vu = \u{a - + ^leU in V, (6) 



2Ngm 



In the Thomas-Fermi approximation, e is small, which 
will be our asymptotic regime. We re-scale the distance 
hy R = d/y/e and define u(r) = Rip{x.) where x = Rr. 
Assume that uj = lu^ and tOy — Xlu with < A < 1 and we 
set D, = fl/eu. Since the trapping potential is stronger 
than the inertial potential, we have 17 < 1/e. The energy 
can be rewritten as: 



E2d{u) 



-l-fj • (m, Vu X r) 



(3) 



Due to the constraint J |itp = 1, we can add to E2D any 
multiple of / jup so that it is equivalent to minimize 

/ |Vup + 2n ■ (m, Vm X r) + -^|it|* - ^a(r)|up 
J ze^ £^ 

where a(r) — a— {x^ + X^y^) for some constant a to be 
determined. Let V be the ellipse {a> 0} — {x^ + X'^y'^ < 
a}. We impose the following constraint on a: 



I «(r) = 1. 
Jv 



(4) 



Indeed, as e tends to 0, the minimizer will satisfy that \u\^ 
will be close to a so that the constraint will be satisfied 
automatically by u if we impose (||). Equation (Q) leads 
to = 2A/7r. If A = 1, that is lo^ — ujy, then I? is a disc 
of radius i?o with Rq = 2/7r. 

To study the problem analytically, it is reasonable to 
minimize the energy over the domain V with zero bound- 
ary data for u. Indeed, when a < 0, the energy is con- 
vex so that the minimizer u goes to zero exponentially 
at infinity (see the numerical observation in Q and the 
analysis on the behavior near the boundary of V as well 
as the decay at infinity of the order parameter in 0). 
Denote H"^{'D) the space of square integrable functions 
defined on the domain V that have square integrable 
derivatives up to the order m, and Hq(T>) the space of 

functions in H^{T>) satisfying the zero boundary condi- 

1/2 

tion. Denote the norm (J!p \v\^') by ||?;|| for any square 
integrable function v, we then consider the problem 



min£'e(w) subject to u G Hq{'D), 



where 



E^{u) = / |Vu|^ + 2n • (m, V' 
Jv 



u X r) 



+ — (a{r)-\u\y. 



(5) 



with u = on dV and /ig is the Lagrange multiplier. The 
specific choice of a in (jj) will imply that the term fi^u is 
negligible in front of au/e^ . 

We want to study the behavior of min E^ [u] as e goes 
to 0. In Section 2, we compute an asymptotic develop- 
ment of the energy and in Section 3, the critical velocities 
of nucleation of vortices and the location of the vortices. 
In Section 4, we study the evolution in imaginary time 
and construct some numerical algorithms. In Section 5, 
we present some computational results and the energy 
diagrams. 



II. ASYMPTOTIC DEVELOPMENT OF THE 
ENERGY 

To study the behavior of the minimizer of the energy 
when e goes to zero, we observe that the form of the en- 
ergy ( P) i s close to the Ginzburg-Landau functional stud- 
ied in]|ll|] , ||l^ where the magnetic field has been replaced 
by a rotating term and similar to |l3| except for the trap- 
ping potential and the minimization over constraint. The 
main idea is to decouple the energy into three terms: a 
part coming from the solution without vortices, a vortex 
contribution and a term due to rotation. The estimate on 
the vortex contribution was developed in 



A. The solution without vortices 

Firstly, we are interested in solutions without vortices, 
that is u has no zero in the interior of T). Thus we con- 
sider functions of the form 77 = Je^^ , where 77 is in i/p (2?) 
and / is real and has no zero in the interior of T). We 
consider first minimizing i?^ over such functions without 
imposing the constraint that the norm is 1, that is, / 
and S minimize 



V 



|V/P + ^(a-/2)2 



j /2|V5 - O X r|2 - fn^r"^. (7) 



We have / = on &D and 



A/ + /V5(V5-2f2xr) = ^/(a-/2)inI?, (8) 



div(/2(V5-Oxr)) = 0. 



(9) 



Equation (|) implies that there exists ^ in H^{'D)r]H^{'D) 
such that 



(10) 
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where V-*-^ = {~dy^, dx£,)- So ^ is the unique solution of 



B. Decoupling the energy 



div( j2 VO = in P, ^ = on dV. (11) 

Note that equation ( |l0|) is the equation for the velocity, 
but we prefer to write it as an orthogonal gradient for 
our later purposes of integration by part. In the special 
case where 2? is a disc, the minimum of (Q) is reached 
for VS = but this is not the case if V is an ellipse and 
there is a non trivial solution of (^. 

a. The case of the disc Assume that A = 1 so that 
I? is a disc and a{r) = i?g — r^. As discussed earlier, 
VS* = in this case so that the energy becomes 

£eif)= ljVf\' + ^iair)-fr. (12) 

Let % be the minimizer of ( p^ over real valued functions 
in Hq{'D). Then, ij^ has no vortex, is independent of fl 
and satisfies 

Arje = -^fleiVe " «) = On dV. 

When e tends to 0, ry^ is close to a except on a boundary 
layer of size e^^^ close to dV. More precisely, using sub 
and super solutions, one can verify that 



v/^tanh((5(v/^)^|log £|) < r/e < 

for \Rq — > Ce^/'^. In fact, one can construct a sub 
solution of the type above in any region |i?Q — | > Ce^/^ 
with P < 2. Then the value of S is less than c(2 — f3). 

The boundary layer can be analyzed using the change 
of variable x = {Rq — r)/e^/^ and Ve{x) — r\^(r) j e^^"^ . 
Then v^, satisfies the Painleve equation 

v" — v{v'^ — 2Rox), v{0) — 0, v{x) ~ 2Rox for x large. 

The boundary behavior has already been studied in ||ic[| 
and 1^, but using matched asymptotics. 

The energy of rye can also be estimated by a test func- 
tion equal to ^/a except on the boundary layer to get 



i?e(%)<y|log£|(l + o(l)). 



(13) 



b. The case of an ellipse As discussed before, the 
minimum rj^ = fgC^^^ of ([7|) has a non-trivial phase, f^ 
tends to a in every compact subset of V and the function 

given by ( |l0|) or (^ij) tends to the unique solution ^ of 



div(-VO = in P, ^ = on dV. (14) 
a 

One can easily get that £,{x,y) = —a'^{x,y)/{2 + 2X'^). 
Using ([l^), we can define Sq, the limit of S^, to be the 
solution of aCV^o-^ri X r) = ilW^^ with zero value at the 
origin. We have 5*0 = Cflxy with C = (A^ - 1)/{X'^ + 1). 
We see that Sq cancels when A = 1 that is in the case of 
the disc. This computation is consistent with the one in 
Pi , though it is derived in a different way. 



Let rye = Ae''^^ be the vortex free minimizer of dis- 
cussed previously without imposing the constraint on the 
norm of u. Let be a minimizer of E,, under the con- 
straint |up — 1 and let Ve = Ui^lri^. Since ry^ satisfies 
the Gross Pitaevskii equation (||)-(p|), we have 



/ {\v\^ - l)(-iA/2 - l/2(a - fl) + |V/,e' 
Jv I e 



~2f^{VSe • O X r) = 0. 

Using this identity, one can get that the energy E^{u^) 
decouples as follows 

E,{u,)^E,{7^,) + G,{v,) 

+2 f \r,,\^{VS, -nxr)- [tv,, Vv,), (15) 

where 

JjVe\'\^Ve\'+^-^{l-\Ve\r- 

This decoupling was used in jl^ in the case of a disc 
where \/Sr = 0. 



C. Estimate on the energy 

We now estimate the terms in (|l5|). The first term 
i?j(77e) is a constant only depending on s, and not on 
the solution type, that is with or without vortices. The 
second term gives a contribution coming from the vortices 
and the third term is due to the vortices and rotation. 

We use the analysis of vortices developed in |l^ and 
later in ||, [l|, |l|. Let V'^ = D \ {x, Aist{x,dV) < 
e^}, with [3 < 1. Then in T>'^ it is possible to define 
vortices for in the following way: there exist balls Bi = 
B{pi,e^ ) where pi are points in T>'^ at mutual distance 
bigger than 8e'' and f3' > f3, such that \vs\ > 1/2 in 
I?[ \ UiBi. Moreover, the degree di — deg{vs/\vs\, dBi) 
is not zero and there is an estimate of the energy of 
in each ball Bi. This analysis means that vortices are in 
fact defined in the balls where is less than 1/2 and has 
a non zero degree. This allows us to compute Ge{ve) for 
which only the gradient term in the vortex balls will give 
a contribution: each vortex gives a contribution in the 
amount of 27r|log e\ due to its degree and a contribution 
of lower order which comes from the interaction with the 
other vortices. Moreover, Iry^p is almost a: 



G,{v) = 27r\\oge\J2 



di\a{pi 



-27r'^didja{pi)a{pj)\og\pi - pj\ +0(1). (16) 



In order to estimate the third term in (15) we let be 
the solution of a(r)(VS'o — $7 x r) = ilW^X^ which is 
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zero on the boundary of V^. That is solves ( p^ but 
with zero boundary data on dV^ instead of dV. Hence 

converges to ^ the solution of (p^. An integration 
by part on the last term of (15) using the definition of 

and the definition of the degree of on vortex balls 
and the fact that the higher order term comes from an 
integration on the vortex balls yields 

/ \ri,\^{'^Se-flxr)-{iv,\/v) 
Jv 

/ Vl-{iv,dX^-KVv){l + o{l)) 

= J2^TTnd,X,{p,){l+0{l)) 

i 

= E T^(" - 1-^1' - >^vri^ + ■ (17) 

i 

Finally, one can derive from (p^-(p^-(p7[) an asymptotic 
development of the energy for a solution with vortices. 

E,{u,) - E,{ri,) ~ 27r|log e\ ^ \dMP^) 

-Pj\ 



-2TT^didja{pi)a{pj)\og\p, 
nil 



1 + A2 



a 



(18) 



Note that the minimal energy for solutions without vor- 
tices in is i<^e(?7e)-|-0(e|log e|): it is not exactly Eg{r]^) 
since % is a minimizer without the constraint ||?7e|| = 1, 
but it almost equals to E^{ri^) since /j, a = 1 and Iry^p 
approaches a asymptotically. 



III. CRITICAL VELOCITIES 

A. Critical velocity for the existence of one vortex 

Let Me be a minimizer of (P) with one vortex at a 
point p inV with coordinates (x, y) and let AE'g be the 
difference between E^iu^) and the energy of a solution 
without vortex {E^{ri^) + 0(£|log ej)): 

(|log e\ ^^{a -x^~ AV)) (1 + 0(1)). (19) 

This expression has been obtained by Svidzinsky and Fet- 
ter 1^ using a different method. The form of Ai?e allows 
the computation of two critical velocities fJs and f2i for 
the existence of vortices: Vis is the velocity for which 
the solution with one vortex starts to be locally stable 
and ill for which it starts to be globally stable. For 
< r^s, Ai?e is a decreasing function of the position 
of the vortex; \p\ = is a local maximum of AiJ^. For 



< ri < fli, \p\ = is a local minimum for Ai?j. Note 
that Ai;e(p e dD) = and /^E^{\p\ = 0) > 0. 

For il > rii, IpI = is the global minimum for AiJ^. : 

- i^llog e\ = ii^x/2^|log e\ , 



2a 



4\/A 
1 + A2 



ill = ii^llog e\ = ^— ^\/2^|log e\ 
a 2VA 



that is 



1 + A2 nh^ , /Ngm\'^n 



-UJ 



4yA V ^9m V n 



o 1 + ^' 

2y/\ 1 



log, 

Ngm V 



( Ngr 



1/2 



Note that Castin-Dum M for the case A = 1 find fli = 

Lu^{TTh^)/{Ngm)\og(^{C/y^){Ngm)/h'^^ , with C ~ 

1.8, hence C/ ^/tt ~ 1 which gives a value of fli very close 
to ours. They also have fii = 2^^ for the case A = 1. 
Increases in anisotropy yield higher fli as already noticed 
in Q, but as A tends to infinity, Qi becomes bigger than 
1/e so that vortices cannot be stabilized. 

It can be proved that there exists which tends to 
zero with e such that for Q < fii — fcg, the minimizer 
of £'e has no vortex and for fl > fti + k^, there exists a 
minimizer with a vortex. Moreover, for fii -I- fc^ < < 
fli + 0{l), any minimizer has one vortex of degree 1 tend- 
ing to the origin. The proof consists in constructing a test 
function with a vortex at the origin and computing the 
energy of this test function. This yields an upper bound 
for the energy. The lower bound relies on estimates for 
Geiv) from and O, ITl. 



B. Critical velocity for n vortices 

Similarly, one can compute , the critical velocity for 
the existence of n vortices. For this purpose, one can 
prove that as e goes to 0, the vortices tend to the origin 
and they are of degree 1 that is rf^ = 1. This is similar 
to ||l^, ijl^. The test function consists in putting the n 
vortices on a polygon centered at the origin of size l/\/f2 
in X and 1 / X^/fi in y. More precisely, we let pi with coor- 
dinates {xi, jji) be such that xi — Xi^/il and yi — Xyi^/fi- 
This allows us to estimate the energy of a solution with 
n vortices centered at pi from (p^ , (p^) , (p7|) : 

Ee{u) = Ee{i]e) + 27rna(|log e\ - ^ .^ ^a) 



+ 7r(n^ — n)Q;^logil -|- w{pi 



1 + A2 
,p„) + a. + o(ll20) 



where C„ is a constant that depends on n and A and 



w(pi, ...,p„) = -7ra^^log(^|xi - Xjf + 



A2 



5 



2 |log e\ 



A2 VLa 



(21) 



Recall that = 2A/7r. For fixed A, w is of order 1, 
hence is of lower order than the previous terms. Then 
the critical velocity for the existence of n vortices can be 
computed from (|20|) 



a„ = (l + A^)J— |l0g£ 



(l + A^) 



(n-l)log((l + A2)J^|log£|),(22) 



and the critical velocity in the original parameters is Vtr, 



(l + A^ , ^ra^ .A^y/. 



2 -"W^^^^H^ 



+ (n-l)i 



-log(- 



27rlog(^)^ 



2NgTrr^\ 2^/A ' ^i' 



C. Location of vortices 

Once f2 is close to f2„, the location of the vor- 
tices is characterized by the configuration of points {pi\ 
which minimizes the function w given by (pl|). In 
non-dimensionalized variables, the points are given by 
_Rpi/^27r|log ej. For convenience, we define 



|log £| 



A V(l + A2) 



Note that given the value of ri„ in (|2^), to leading order, 
p is equal to \/27r/A/ (1 + A^). We use the value of a and 
p to get a simplified expression for w. 



w{pi, ...,p„) = -2\{^\og(^X^ - Xj\^ + ^2 



(23) 



The critical points of w, and thus the vortex positions, 
satisfy: 



E'^"'(5i ^ S?) 
■ 



A2|ij - + \y^ _ ^j-p 
- 



An immediate observation is that 



(24) 
(25) 

(26) 



By multiplying the equations with Xi and iji respectively 
and adding the results together, one can obtain 



J2iS^! + yf) = nin-l)/i2p) 



(27) 



Similarly, multiplying the equations with and — A^Xj 
respectively and adding the results together, one gets 



p(l-A2)^i,y, = 0. 



(28) 



Unlike equations ( [2^p7|) where the dependence on A is 
implicit, the equation (|28|) leads to a property 



E x.,yi = , for A 7^ 1 . 



(29) 



The above observations lead to more precise predictions 
on the location of vortices. 

For instance, in the case n = 2 we get ii = —X2 and 
yi = — j/2- For A = 1, we have an infinite set of solutions 
consisting in two points on the circle xf + yf = 1/p , 
symmetric with respect to the origin. For A 7^ 1, ( [29| ) 
leads to Xijji — for i — 1,2, and we have a pair of 
solutions with either Xi = 0, jji = ±a/1/2p, or jji = 0, 
Xi = ±^yi/2p. Checking the corresponding values of w, 
we get that for A 7^ 1 , the minimizer of w corresponds to 
having both vortices staying on the long-axis of the ellipse 
in the original scaling (that is on the x axis if A > 1 and 
on the y axis if A < 1). The estimate of the location is 
in agreement with the numerical solutions given later. 

For the case n = 3, we also get that the 3 vortices are 
on the long-axis of the ellipse: one centered at the origin 
whereas the other pair stays symmetrically on the long- 
axis with Xi = zty^3/2p if A > 1. Similar discussions can 
be carried out for other values of n. In general, for the 
case of n vortices, we let — n{n—l)/2p and Xi — RiXi, 
iji = Riyi- It follows from ( p7| ) that the points {pi} are 
localized by the minimum of 



Elog(|. 



\Vi~yj\ 

A2 



under the constraint Y^iXj + yf = I. It is reasonable to 
expect that this formula will lead to a vortex array as ob- 
served in p^ . Note that the minimizer {pi = RiPi} has 
no explicit dependence on O, thus, we expect for a given 
A, the vortex configuration to be of similar structure for 
values of f2 close to ri„. On the other hand, for a given n 
and ri, we expect, however, that for sufficiently large A, 
that is for highly anisotropic traps, the minimizer of w is 
given by the coUinear solutions with vortices all located 
on the long axis of the elliptical trap. 

To our knowledge, theoretical investigation had been 
restricted to the critical velocity for nucleation of one 
vortex. Here, we are able to deal with the case of multi- 
ple vortices and to precisely characterize the location of 
the vortices. Our estimates are also consistent with the 
numerical results given later. 
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IV. EVOLUTION EQUATION AND 
NUMERICAL SCHEMES 

To numerically compute the energy minimizers of (|^), 
we notice that the energy in (o) can be rewritten as 



\{y-iA)u\^ + —{aM~\u??\+c, (30) 



where ae(r) = a(r) — e^il^r^ 



f2, and c, = Jj^{a^{r)^al{r))}. 



The above formulation of the energy has a striking sim- 
ilarity with the high-kappa high-field Ginzburg-Landau 
energy with a variable coefficient p^ . 



A. Evolution in the imaginary time 

To numerically compute the minimizers of (|^) , we con- 
sider the time-dependent equation in the imaginary time: 



du 
'dt 



(V-iA) u- 



|^^P^^- 



a£(I■), 



^J,eiu)u (31) 



in V with initial condition M(r,0) = uo(r) in V and 
boundary condition u = on dV. Here, fieiu) denotes 
the Lagrange multiplier. Assume that uq satisfies the 
constraint ||mo|| = 1- Then, by taking 



we get 



1^ ||(v - ,:a) + - ^H^l , 



-1=0 



Thus, the constraint jitp = 1 is ensured at all time. 
Moreover, using {u,Ut) = 0, we get the energy estimate: 

Thus, we easily get that for any (0,T), if uo G Hq{T>) 
and ||uo|l = 1, there exists a unique strong solution u of 
(^) satisfying the constraint = 1. Using argument 
similar to that in [ pT| , we may also get that ast ^ oo, 
u approaches to a steady state solution which is a critical 
point of the energy. Well-posedness for initial data 
may also be be obtained. 



high-kappa high-field time-dependent Ginzburg-Landau 
equations (20) , and adapt a code developed in jTs], |l^, |o) . 
Spatially, we use a standard finite element approxima- 
tion, see [|l^, |l^ for details. Here, we focus on the time- 
discretization and the treatment of constraint. It has 
been observed that there are some steady states exhibit- 
ing meta-stability, thus it is important to get asymptot- 
ically stable schemes for large time which in general re- 
quiring the use of implicit schemes with no limitations 
on the time step size. 

Let {u„} be approximate solutions of {u{tn)} at dis- 
crete time {tn} with time-step At„ = t„ — tn-i- We 
discuss two time-discretization schemes and also some 
results of numerical experiments. 

c. A first order backward- Euler in time discretization 
Given Un-i, we first solve for u*: 



— ( V — lA) U — fl(Un^i)U 



Atn 



(32) 



Then, we apply the projection u„ = Both the 

backward Euler step and the projection step gives only 
first order in time accuracy. 

d. A norm-preserving, energy-decreasing second or- 
der scheme For any u, v and their complex conjugate 
M, V, we let f{u,v) — -I- |wp)(u + t')/2 which satis- 
fies f{u, v){u — v) — (\u\* — |u|'')/2. Given m„-i, we first 
solve for u*: 

2[U* — Un-l) . . n2 * / *N * 

^ — — [y — iK) u — v{u )u 

-|-^/(2m* - w„_i,u„_i) - ^fleU* = (33) 
where v{u*^ is given by 

Ku*)^KI' = ^{l(v-zA)«*|2} 

+ / I — /(2u* - u„_i, u„_i)t2* - 



Then we let u„ = 2u* — Un-i- Taking the inner product 
of the equation (|3^ ) with u*, we get (u* — u„-i, u*) = 0, 
which leads to ||u„|p = ||Mn-i|P- That is, the norm is 
preserved at each time step. Taking the inner product of 
the equation (p3|) with u* — u„_i, it is easy to get 



-I- £e(w„) - Ee[Un-\) = 0. 



B. Numerical schemes 



There are various ways to solve the time-dependent 
Gross-Pitaevskii equations, see for example or Q. We 
take the advantage of the striking similarity with the 



Thus, during the discrete time evolution, the energy de- 
creases. This discrete scheme is second order in time 
and unconditionally stable. It also preserves some essen- 
tial features of the continuous dynamic system, making 
it suitable for long time integration and for studies of 
meta-stabilities of the solutions. 
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C. Description of the numerical experiments 

We have used the above schemes to calculate various 
numerical solutions for the parameter values e — 0.02, 
A = 1 and A = 1.5. The spatial finite element space is 
taken to be C° piecewise quadratic elements on trian- 
gular meshes. As we are mostly interested in the mini- 
mizers of the Gross-Pitaevskii energy, the time-evolution 
is employed as a means of marching to the steady state 
solutions. For this reason, we have used variable time- 
steps in order to accelerate the convergence in time. The 
nonlinear systems are solved by a Newton like methods 
at each step. Though such a method is computational 
costly per step, this drawback is offset by its uncondi- 
tional asymptotic stability for marching to the steady 
state. We have also computed the solutions using refined 
meshes to ascertain the numerical convergence. 

To obtain solutions for various velocities, we have used 
a number of differential initial conditions. For example, 
we have used |uo(r)P — a(r) for r G I? which serves 
as a good approximation to the steady state solution, 
especially in the case of vortex- free solutions. We note 
that for large values of 51, this choice of initial condi- 
tion can also lead to steady state solutions with multiple 
vortices. Detailed solution branches are described in the 
next section. In addition, we have also used other initial 
conditions that manually plant vortices in the domain in 
order to find different solution branches. Finally, a con- 
tinuation in the parameter fl has often been employed to 
follow a particular solution branch and to compute the 
bifurcation diagrams. The continuation procedure also 
provides a test for the local stability of the numerical so- 
lutions: when one branch becomes unstable, the solution 
jumps onto a different branch. 

V. NUMERICAL RESULTS AND 
BIFURCATION DIAGRAMS 

We now present some pictures of numerical solutions 
and discuss the various solution branches. 



A. Description of solutions 

e. The case of a disc For any f2, there is a vortex 
free solution, which is close to a{x, y) except near the 
boundary. For = 0, in addition to this vortex free 
solution, there is also a one vortex solution, as illustrated 
in the first column of Figure ^. 

For larger fl, solutions with multiple vortices are shown 
in the other columns of Figure |l| For instance, solutions 
with 2 and 3 vortices for U, — Ih are shown in the second 
column, solutions with 4 and 5 vortices for Q. — 17.5 are 
shown in the third column, and solutions with 6 and 7 
vortices for 17 = 20 are shown in the last column. The 
parameter values for which the single-vortex and multi- 
vortices solutions exist are to be presented in the next 




FIG. 1: Contour plots of \u\ for = (1st column), 15 (2nd 
column) 17.5 (3rd column) and 20 (4th column). 



section. Figures ^ and [3| provide the surface plots (and 
better view) of |m| for a vortex free solution at SI = 
and a solution at f2 = 20 with four vortices respectively. 
Each solution has a top and bottom view, the paraboloid 
shapes are easy to visualize from the picture. 




FIG. 2: Surface plots of \u\ for solutions at f2 = 0. 

/. The case of an ellipse We now present some so- 
lutions for an ellipse corresponding to A = 1.5. In Figure 
^, the contour plots of the magnitude of the solution |it| 
with ri = 17.5 are drawn while in Figure |, S7 = 25. 
For S7 — 17.5, it is interesting to compare the location 
of vortices with the analysis of section III.c: in this case 
p = 0.79, the long axis is 0.9887 and we find for the lo- 
cation of the vortex x = 0.19 for n — 2 and x — 0.33 
for n = 3. The picture shows that the vortex is about 
1/5 of the long axis for 2 vortices and 1/3 for 3 vortices, 
which is consistent with the analysis. We believe that for 
a bigger number of vortices n, their location corresponds 
to the positions minimizing (P3[). 

Note that there were two vortex configurations with 
five vortices for D. = 25, among which, the non-symmetric 
configuration corresponds to a solution with a lower free 
energy, though the difference in their energy values is 
very small. 
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FIG. 3: Surface plots of \u\ for solutions at 57 = 20. 




FIG. 4: Contour plots of \u\ for Q = 17.5, A = 1.5 



B. Branches of solutions 

An issue that we have studied is the existence of 
branches of n vortex solutions as Q is varied and espe- 
cially which one is the minimizer. 

g. The case of the disc For 17 = 0, the solution with 
lowest energy is the vortex free solution. We start with 
this solution as initial value for the time dependent prob- 
lem for a slightly bigger fi. This device allows us to 
continue the branch of the vortex free solution as is in- 
creased. We find that the vortex free solution is obtained 




FIG. 5: Contour plots of \u\ for = 25, A = 1.5 



as the limit when t is large of the evolution equation up to 
n = 19. For = 19, six vortices are nucleated from the 
boundary and eventually one vortex moves to the center 
and final configuration is similar to that in Figure ^ (top 
right). Now if fl is decreased from the value 19 using 
the 6 vortices solutions as initial value, we see that this 
6 vortices solution branch exists down to O = 16 when 
it jumps to 4 vortices. If we decrease fi further, then 
we stay on the 4 vortex branch down to il = 13 when it 
drops to a 2 vortex branch. Similarly, if we increase f2 
from 13, the 2 vortices solution will persist up to = 21. 

As for the one-vortex solution branch, it is computed 
by planting a vortex-like function at the center in the 
initial condition for il = 10, then the branch is computed 
by continuation in (both increasing and decreasing). It 
is interesting to observe the fact that this branch extends 
all the way to H = 0, such persistence of the one- vortex 
solution even for the zero velocity has been elaborated 
by various authors, see Q| for instance. Since implicit 
integration is used, we are able to indeed march to to the 
steady state and to ascertain that this persistence is not 
due to the metastability. On the other hand if the vortex 
is planted away from the origin, it disappears for small 

as we will see later. 

h. The case of an ellipse The same kind of be- 
haviour occurs for the vortex free solution branch with 
A = 1.5: when Q is increased from 0, we stay on the vor- 
tex free branch up to 17 = 22.5 when the solution jumps 
to the 4 vortex branch of Figure ^j. If we decrease fl 
starting from this 4 vortex branch, the solution will stay 
on it down to 17 = 15 when it jumps to 2 vortices and at 
17 = 10 it jumps to the vortex free solution. Similarly to 
the case of the disc, if one vortex is planted at the center 
at time 0, it will persist in time even down to 17 = 0. 

For 17 large enough, several vortices are nucleated at 
the same time. Basically, we are counting on initial con- 
ditions and the round-off errors to break the symmetry 
of the region and strong symmetry presence often makes 
symmetry breaking much harder to achieve. For the disc 
case, we expect equal chance of vortex nucleation from 
any point of the boundary. It turns out that for 17 — 19, 
an unstable front produces oscillations of almost equal 
magnitude, and spins off 6 vortices at the same time. 
Had the disturbance being unevenly distributed, then it 
would be possible for some vortices to get spin off ahead 
of others, and due to the repulsion, the others may never 
have a chance to appear, thus we may see solutions with 
a smaller number of vortices. 



C. Energy diagrams 

We now discuss the energy diagrams in relation to the 
discussion of the critical velocities given in the earlier 
sections. In Figures ^ and we have plotted the energy 
given by (0) as a function of 17 for the various branches of 
solutions (according to the number of vortices). Again, 
e = 0.02 in our computation. 
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FIG. 6: The energy vs. Q curves for A — 1.0 



'2vortices' 
'3vortices' 
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FIG. 7: The energy vs. Q curves for A — 1.5 



i. The case of the disc (Figure]^) As discussed ear- 
lier, the vortex-less solution, in the case of the disc, ex- 
ists for all values of f2 and is independent of fi, we thus 
have a constant line for its energy. These vortex-free 
solutions are the global energy minimizers for small f2 
{ft < 9.3) whereas for Q > 9.3, they have larger energy 
than the one- vortex solution. For multi- vortex solutions, 
we see that each becomes the global energy minimizer for 
a range of values of Q. It is interesting to compare this 
result with the value found in section 2 where ili — 9.8. 
Similarly, we obtain from ( p^ ) that fin — fln^i — 1.77. 
Though the value of fli is slightly overestimated, the dif- 
ference fin — fln-i looks good for small n: the numer- 
ics indicate fl2 = 12.0, = 13.6, fl4 = 15.8 and our 
theoretical computations yield ft2 — 11.6, fl^ = 13.4, 
fls = 15.2. 

However, when fl is increased from 0, we saw that we 
stay on the vortex free branch up to 57 = 19. This means 
that the vortex free solution is a local minimizer up to 
this value. We do not have any theoretical estimate for 
this value of local minimum. Similarly, for multi-vortex 
solutions, hysteresis loops are present. For the solution 
with six vortices, there are two possible configurations, 
one with all six on a concentric circle, one with only five 
on a concentric circle while the remaining vortex at the 
center of the disc. This occurs for fl — 22.5. Though the 
difference in the energy is hardly noticeable, the solution 
with a center vortex does have a smaller energy value. 

j. The case of an ellipse (Figure^ Now, we discuss 
the elliptical case with A = 1.5. The energy versus fl 
curves are given in Figure |^. 

The vortex-less solution, in the case of the ellipse, is 
no longer independent of fl as in the case of the disc, 
as illustrated by the dependence its energy on Q. Aside 
from that, the hysteresis phenomena also occurs just like 
for the case of a disc. 

Our numerics indicate that Qi = 12, fl2 ^ 15. It is 
interesting to compare this result with the value found 
in section 2: we obtain that fli ~ 13 from ( |2^ ) that 
fin — fln-i — 2.88. Though the value of fli is slightly 
overestimated, the difference fln — fln-i looks good again 



since we find fl2 — 14.9. 



D. Displacement of the vortex from the center 

Based on the earlier estimate (^9|) for the one-vortex 
solution, we see that for small fl, the displacement of the 
vortex away from the center leads to the drop of energy 
{fl < 4.9). For slightly bigger SI a vortex at the center 
is a local minimum (fl < 9.8). Let us analyze the time 
dependent equation using a vortex off center as initial 
condition and let us examine how it evolves. We find 
that in the case of the disk, for fl < 7, a displacement 
on the size of a tenth of the radius causes the displaced 
center vortex in the one- vortex solution to move towards 
the boundary. For fl > 10, the center vortex in the one- 
vortex solution moves back towards the center if under 
a displacement on the size of a third of the radius. For 
intermediate values of fl, the vortex moves back to the 
center for small displacement but moves towards bound- 
ary for large displacement. 

The following two pictures (Figures ^ and ^ show, in 
the case of a disc and a small displacement, the marching 
away of the center vortex for the solution with f2 = 
(starting from top row then to the bottom row) and the 
marching back to the center for the solution with fl =^ 7.9. 
















FIG. 8: Perturbing vortex away from center for S7 = 
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FIG. 9: Perturbing vortex away from center for Q = 7.9 

VI. CONCLUSION 

We have presented a new framework for the study of 
the Gross Pitaevskii energy in the Thomas Fermi Hmit: 
we have defined a nondimensionnahzcd parameter e and 
we have computed theoretically an asymptotic develop- 
ment of the energy, the critical thermodynamic veloc- 
ities of nucleation of vortices and the location of vor- 
tices. This extends the results of We have proposed 
and implemented time integration schemes which enjoy 
the norm-preserving and energy-descreasing features and 
thus useful for studying the stability and metastability of 
solutions. We have also presented energy diagrams com- 
puted numerically for the various vortex solutions. We 
have noticed that our theoretical predictions for critical 



thermodynamic velocities are quite consistent with the 
numerics which encourages us to think that our approx- 
imation £ small is correct though |log e| is not small. 

In our computation, we have taken e = 0.02 and the 
number of vortices we have observed ranges from to 10 
for fl between and 25. 

Our aim for a future work is to use the results of 
for Ginzburg-Landau in dimension 3 to get results for 
rotating Bose Einstein condensates in dimension 3, and 
especially the critical velocities and an asymptotic devel- 
opment of the energy. In particular, it is obtained in 
that Qg = 3/5fli for a vortex line parallel to the z axis, 
but we would like to get an expression of the energy for 
any type of vortex line. We will also carry out further 
numerical simulations in 3 dimension to compare with 
experimental data. 



Acknowledgments 

The authors would like to thank Y.Castin for very 
interesting discussions and also T. Riviere, E. Sandier, 
S.Serfaty. This work is supported in part by a joint 
France-Hong Kong research grant. 



D.Butts and D.Rokhsar, Nature 397, 327 (1999). 
Y.Castin and R.Dum, Eur. Phys. J. D, 7, 399 (1999). 
F.Dalfovo, S.Giorgini, L. Pitaevskii and S.Stringari, Rev. 
Mod. Phys. 71,463 (1999). 

D.L.Feder, C.W.Clark and B.I.Schneider, Phys. Rev. 
Lett., 82, 4956 (1999). 
A.L.Fetter, Phys. Rev A, 148, 429 (1965). 



A.L. Fetter and A.A.Svidzinsky, |cond-mat/010200; . 

K. Madison, F. Chevy, W. Wohlleben and J. Dalibard, 

Phys. Rev. Lett., 84, 806 (2000). 

A.A.Svidzinsky and A.L.Fetter, Phys. Rev. Lett., 84, 
5919 (2000). 

F.Dalfovo, L. Pitaevskii and S.Stringari, Phys. Rev. A 54, 
4213 (1996). 

A.L.Fetter and D.L.Feder, Phys. Rev. A, 58, 3185 (1998). 
A.Aftalion, E. Sandier and S.Serfaty, to appear in J. 
Math. Pures et Appl. (2000). 

S. Serfaty, Comm. Contemporary Mathematics, 1, 213 
and 295 (1999). 

S.Serfaty, to appear in Control, Optimization and Calcu- 
lus of Variations (2001). 



[14] F. Bethuel, H. Brezis and F. Helein, Ginzburg-Landau 

Vortices, (1994) Birkhauser. 
[15] F. Bethuel and T. Riviere, Annales IHP, Analyse non 

lineaire, 12, 243 (1995). 
[16] D.L.Feder, C.W.Clark and B.I.Schneider, Phys. Rev. A, 

61 011601(R) (1999). 
[17] S. Chapman, Q. Du, M. Gunzburger and J. Peterson, 

Adv. Math.Sci. Appl. 5, 193 (1995). 
[18] Q. Du, M. Gunzburger and J. Peterson, Phys. Rev. B, 

46, 9027 (1992); Q. Du, M. Gunzburger and J. Peterson, 

Phys. Rev. B, 51, 16194 (1995). 
[19] Q. Du, M. Gunzburger and J. Peterson, SIAM Review, 

34, 54 (1992). 

[20] Q. Du, P. Gray, SIAM Appl Math, 56, 1060 (1996). 

[21] F. Lin, Q. Du, SIAM Math Anal, 28, 1265 (1997). 

[22] L. Simon, Annals of Math 118, 525 (1983). 

[23] F.Pacard and T. Riviere, in Progress in Nonlinear Differ- 
ential Equations and their Applications, 39 Birkhauser 
Boston, Inc., Boston, MA, 2000. 



